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ABSTRACT 


In  this  work,  several  algorithms  based  on  higher-order  moment  (autocorrelation) 
matching  of  single  hydrophone  element  data  have  been  developed  and  tested  on  real 
transient  data  sets.  Of  particular  interest  is  the  success  and  robustness  of  the  Frequency- 
domain  Autocorrelation  Matching  (FACM)  algorithms  in  the  presence  of  environmental 
mismatch,  signal  mismatch,  and  noise,  for  different  signals  in  an  unknown  environment. 
Recently  acquired  data  was  analyzed  for  signal  variability  in  terms  of  spatial  coherence 
of  phones,  beams,  and  modal  structure.  The  ability  to  localize  using  these  higher-order 
moment  matching  algorithms  was  compared  to  the  spatial  structure  of  the  signal,  the 
placement  of  the  receiving  elements,  and  the  signal  variability. 

This  work  suggests  that  the  FACM  algorithms  are  strongly  dependent  on  the 
source-receiver  relative  positions,  and  on  the  uniqueness  of  the  signal  vertical  structure. 
It  is  also  shown  that  their  performance  increases  with  the  number  of  multipath  arrivals 
and,  therefore,  with  the  range.  More  importantly,  the  localization  results  obtained  with 
raw  linear  frequency  modulated  (LFM)  signals  seemed  to  be  as  useful  as  the  ones 


obtained  from  matched-filtered  data. 
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I.  INTRODUCTION 


In  the  past,  passive  sonar  systems  were  only  able  to  provide  bearing  information,  and 
the  estimation  of  source  range  relied  on  analytical,  over-simplified  techniques.  Starting  in 
the  mid-60’s,  however,  with  relatively  advanced  computational  resources,  the  received 
signals  have  been  processed  in  order  to  also  provide  information  about  the  distance  to  the 
source.  Today,  one  of  the  most  used  and  successful  research  techniques  is  matched-field 
processing  (MFP)  where  a  comparison  is  made  between  the  parameters  6f  the  received 
signal  and  those  of  signals  generated  by  a  synthetic  source  virtually  positioned  at  each  point 
in  the  search  grid.  A  match  is  obtained  when  the  correlation  of  the  synthetic  and  the  true 
signals  are  at  a  maximum. 

Transient  signals  are  especially  important  for  naval  purposes,  since  they  can  represent 
several  instances  of  a  ship’s  routine,  such  as  firing  a  torpedo,  opening  of  hatches,  starting 
of  pumps,  etc.  Therefore,  the  Transient  Localization  Project  at  the  Naval  Postgraduate 
School  has  been  studying,  since  1993,  several  algorithms  for  localizing  transient  sources 
based  on  methods  similar  to  MFP. ' 

In  1995,  three  distinct  schemes  were  studied:  (a)  the  fully  coherent  localization,  which 
was  considered  impractical  as  it  required  an  extremely  precise  (and  unavailable) 
propagation  model,  along  with  accurate  environmental  descriptions;  (b)  the  semi-coherent 
localization,  which  correlated  only  the  slower  varying  amplitude  terms  of  the  signal, 
trying  to  surpass  the  limitations  of  method  (a);  and  (c)  the  autocorrelation  matching,  which 
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compared  the  autocorrelations  of  the  received  and  modeled  signals  using  a  simple  Bartlett 
technique.  The  NPS  results  showed  that  method  (b)  was  not  successful  in  localizing  the 
signals  of  naval  interest,  and  method  (c),  while  not  perfectly  consistent,  was  the  most 
useful  technique.  According  to  Miller  et  al.  (1996),  autocorrelation  matching  has  an 
advantage  over  the  other  methods,  since  its  processing  allows  the  notching  of  the  noise 
energy  and  the  possibility  of  a  multi-transient  gain. 

A.  THE  NUWC-TSP  EXPERIMENT 

The  transient  signal  processing  (TSP)  experiment  was  a  field  study  sponsored  by 
NUWC,  and  conducted  both  by  NUWC  and  C&M  Tech.,  Inc.  on  28-31  July  1997.  Its 
purpose  was  to  provide  data  to  study  the  influence  of  a  highly  variable  shallow  water 
environment  on  acoustic  propagation,  and  the  effectiveness  of  localization  algorithms. 

The  study  area  is  shown  in  Figure  1.1  and  as  it  is  influenced  by  oceanographic  and 
geological  processes  related  to  the  shelfbreak,  it  provides  a  complex  coastal  environment 
for  sound  propagation.  The  continental  shelf  water  is  cold,  and  it  presents  large  variations 
with  location  and  season.  The  slope  water,  on  the  other  hand,  is  warmer  and  influenced  by 
the  Gulf  Stream. (Miller,  C.  W.,  1998) 
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Synthetic  signals  were  generated  and  transmitted  in  a  variety  of  scenarios  which  were 
differentiated  by  the  source  depth,  source-receiver  range,  and  bearing.  Two  projectors 
(ITC-104  and  HX-188),  a  32-element  vertical-line  hydrophone  array  (VLA),  and  data- 
recording  equipment  constituted  the  main  assets  of  the  experiment.  The  YLA 
configuration  is  depicted  in  Figure  1.2,  and  the  experiment  configuration  for  one  of  the 
scenarios  is  shown  in  Figure  1.3. 
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VLA  configuration 
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Figure  1.2.  The  VLA  Configuration. 


Figure  1.3.  The  NUWC-TSP  Experiment  Configuration. 
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It  is  important  to  note  that  hydrophones  17  and  18  were  non-operative  throughout  the 
experiment.  In  this  work,  we  were  only  interested  in  the  linear  frequency  modulated 
(LFM),  short  duration  (100  ms)  signals  -  the  transmitted  chirp  and  a  sample  of  a  received 
burst  can  be  seen  in  Figure  1.4.  Their  frequency  spectra  are  observed  in  Figure  1.5. 

The  processed  data  corresponded  to  three  particular  cases:  scenario  9  (source  located 
at  4650  m  range,  45.7  m  depth),  scenario  11  (source  at  460  m  range,  45.7  m  depth),  and 
scenario  16  (source  at  2194  m  range,  27.4  m  depth).  The  information  on  source  location 
was  provided  by  NUWC,  based  on  ship’s  log  and  GPS  data.  While  some  uncertainties  are 
present  in  this  information,  inherited  from  the  processes  used  to  gather  it,  we  assumed  it  as 
our  reference.  The  recorded  tapes  also  contain  array  tilt  information,  which  was  neglected 
in  this  work  for  simplicity. 


Figure  1.4.  The  Transmitted  LFM  Burst  (Left)  and  the  Received  Signal  at  VLA  (Right)  in 
the  Time-Domain.  Note  That  the  Axis  Are  Different  on  Each  Plot 
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Figure  1.5.  Frequency  Spectra  of  Transmitted  Chirp  (Left)  and  Received  Burst  (Right). 

In  order  to  support  the  analysis  of  the  experimental  data,  six  conductivity-temperature- 
depth  (CTD)  casts  were  made.  Two  of  them  did  not  present  useful  data  and  were 
disregarded.  The  other  four  provided  the  sound  speed  profiles  (SSP)  shown  in  Figure  1.6, 
which  also  displays  the  date  and  time  (Zulu)  they  were  taken.  For  processing,  however,  we 
used  only  SVj97_08  (applied  to  scenario  9),  and  SVj97_12  (applied  to  scenarios  11  and 
16),  because  they  were  gathered  at  times  closer  to  those  of  the  studied  runs.  Still,  it  is 
important  to  note  the  large  variations,  particularly  near  the  lower  part  of  the  thermocline 
between  20-35  m.  The  associated  pycnocline  supports  the  propagation  of  nonlinear  soliton 
waves  which  travel  from  the  edge  of  the  slope  onto  the  shelf.  Such  oceanographic  features 
can  create  localized  regions  of  high  sound  speed  contrast  as  strong  as  25  m/s  over  a  100  m 
horizontal  range  in  the  direction  of  the  soliton  propagation  (Chiu  et  al.,  1997). 


6 


Figure  1 .6.  Sound  Speed  Profiles  (SSP’s)  Collected  During  Experiment. 

B.  THESIS  OBJECTIVE  APPROACH  AND  OUTLINE 

The  objective  of  this  thesis  is  to  examine  the  influence  of  shallow  water  acoustic 
variability  on  frequency-domain  autocorrelation  matching  (FACM)  localization.  In  order 
to  perform  such  analysis,  we  processed  the  acoustic  VLA  data,  identified  the  LFM  bursts, 
and  high-pass  filtered  them  to  reduce  the  array  strum  influence  (cut-off  at  50  Hz).  The 
resulting  bursts,  as  well  as  the  replicas  generated  by  the  MMPE  propagation  model,  were 
applied  to  a  Bartlett-type  processor  based  on  matching  broadband  signal  autocorrelation 
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localization  algorithm).  The  acoustics  variability  study  followed  with  the  plane-wave 
beamforming,  modal  amplitude  structure,  and  spatial  coherence  analyses. 

The  remainder  of  this  thesis  consists  of  five  chapters.  Chapter  II  describes  the 
Monterey-Miami  Parabolic  Equation  (MMPE)  propagation  model.  Chapter  III  discusses 
the  Frequency-domain  Autocorrelation  Matching  (FACM)  localization  algorithm.  Chapter 
IV  describes  the  techniques  used  to  evaluate  the  environmental  influences  and  process  the 
signals.  Chapter  V  presents  the  localization  results.  Chapter  VI  concludes  this  study  with  a 
summary  of  findings  and  recommendations  for  future  work. 


II.  THE  MONTEREY -MIAMI  PARABOLIC  EQUATION 
PROPAGATION  MODEL 


The  Monterey-Miami  Parabolic  equation  (MMPE)  acoustic  propagation  model 
(Smith,  K.,  1996)  was  used  to  generate  the  replicas  which  simulated  the  received  signals 
from  a  synthetic  source  on  the  gridpoints  of  a  2-D  search  space  (depth  x  range).  It  allowed 
us  to  predict  the  arrival  structure  at  the  VLA  due  to  a  transient-like,  broadband,  point 
source.  In  this  chapter,  we  introduce  the  general  theory  behind  the  parabolic  equation 
model,  as  well  as  the  method  used  for  its  implementation  -  the  Split-Step  Fourier  (SSF) 
method  (Hardin  and  Tappert,  1973). 

The  parabolic  equation  method  for  computing  underwater  acoustic  propagation  was 
introduced  by  Tappert  (1977).  To  begin  the  parabolic  approximation  development, 
consider  the  representation  of  the  time-harmonic  acoustic  field  in  cylindrical  coordinates 

P(r,z,ty,'f,t)  =  pj(r,z,ty)e~l2nft .  (1) 

As  the  model  is  based  on  an  approximation  to  the  Helmholtz  wave  equation,  substituting 
Eq.  1  into  the  latter  in  .cylindrical  coordinates  leads  to  (Jensen  et  al.,  1994) 


where  n(r,z,  <b)=  — ~ — r  is  the  acoustic  index  of  refraction,  c„  is  the  reference  sound 

c(r,z,  <j>) 

speed,  c(r,  z,  ( |>)  is  the  acoustic  sound  speed,  and  k0-  —  is  the  reference  wave  number. 

Co 

In  this  derivation,  however,  we  neglected  density  variations,  which  could  be 
incorporated  into  a  new  index  of  refraction  without  any  loss  of  generality.  Note  also  that 
the  environmental  characteristics  are  within  c(r,  z,  (j>) ,  and  the  reference  pressure  level 
P0  is  defined  as  the  pressure  amplitude  at  a  distance  R0  -  1 . 

By  assuming  the  ocean  acts  as  a  waveguide  with  a  cylindrical  coordinate  system, 
acoustic  energy  is  primarily  propagated  outward  from  a  source  in  the  horizontal  direction. 
The  pressure  field,  therefore,  can  be  approximated  by 

p/j,Z,  cj>)  =  Vf(r,z,.ty)HW(k-0r),  (3) 

where  the  slowly  varying  envelope  y^(r,z,  ( j>)  modulates  the  outgoing  zero-th  order 
Hankel  function  of  first  kind  H^\k0r). 

In  the  far  field  ( kQr  »  1 ) ,  we  can  approximate  the  Hankel  function  by  the  first  term  of 
its  asymptotic  expansion  (Gradshteyn  and  Ryzhik,  1994), 


H£Xk0r)  = 


(4) 
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and  Eq.  3  becomes 


(5) 

Note  that  Eq.  5  is  normalized  such  that  at  r  =  R0 ,  |\|t|  =  1  and  |py|  =  PQ.  It  expresses 

the  relationship  between  the  PE  field  function  z,  <())  and  the  acoustic  pressure 
pfj,z,  <j)).  By  substituting  it  into  Eq.  2  and  dropping  the  source  term  on  the  right  hand 
side,  we  obtain 


Pf(r,z,i !>)  =  P0  —  \|/(r,z,<}>). 

W  T  f 


If  we  neglect  the  azimuthal  coupling  and  the  far  field  terms,  Eq.  6  can  be  written  as 


r?.  =  Vc+U'p)'*' 


k0Br 


<7) 


where 


T 

op  2  JL 


WzV 


(B) 


and 


(9) 


Note  that  this  is  only  valid  for  once  we  used  the  uncoupled  azimuth 

o(()2 


approximation. 
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Eq.  7  is  known  as  the  standard  parabolic  equation  (SPE),  and  the  accurate  solutions  are 
limited  to  a  half  beam  width  of  10°  to  20°  for  the  propagation  angle.  However,  in  order  to 
extend  this  limit  to  40°  -  70°  ,  a  higher  order  wide-angle  parabolic  equation  (WAPE) 
approximation  (Thomson  and  Chapman,  1983)  can  be  used.  It  also  features  less  sensitivity 
to  the  choice  of  k0 ,  and  to  the  phase  errors  in  typical  deep  ocean  conditions,  relative  to  the 

SPE  approximation.  (Jensen  et  al.,  1994  and  Chin-Bing  et  al.,  1993) 

The  WAPE  operators  are  defined  as 


twape 


1  a2 

k2dz2 


1  + 


2\ 


1JL 

k2dz2 


(10) 


and 

Uwape  s  “ 1 (n  —  1 )  •  (11) 

In  order  to  numerically  solve  the  parabolic  equation,  the  MMPE  uses  the  split-step  Fourier 
(SSF)  method.  This  algorithm  integrates  the  solution  in  range  by  applying  the  operator 
U wape  *n  the  z-domain.,  and  the  operator  TWAPE  in  the  kz  -domain.  The  latter  is  defined, 

in  the  kz  -domain,  as 


TwAPE(k )  =  1  - 


kD 

i--S 

KoJ 


(12) 


Note  that  both  operators  are  just  scalar  multipliers,  and  may  be  applied  independently. 
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The  approximate  solution  for  \| /  is  iterated  in  range  according  to  the  following 
expression 

V(r  +  Ar,2)  =  (13) 

where  the  FFT  calculation  is  consistent  with  the  form  described  in  (Press  et  al.,  1988).  A 
forward  FFT  and  double  conjugation  emulate  the  inverse  FFT.  The  model  outputs  the  field 
functions  \|/  (both  magnitude  and  phase),  referenced  to  a  unit  magnitude  at  r  =  lm,  and 
computed  at  the  spatial  grid  points. 

The  source  assumed  in  our  model  had  a  600  Hz  bandwidth  centered  at  900  Hz. 
Therefore,  the  acoustic  field  was  computed  for  each  one  of  the  discrete  frequencies 
(N=512)  of  the  considered  BW,  and  the  model  properly  represented  the  broadband 
acoustic  field. 

As  the  single  component  of  the  general  field  can  be  expressed  as  Eq.  5,  and  assuming 
a  windowed  (Hanning),  normalized  source  amplitude,  the  time-domain  complex  pressure 
field  can  be  written  as 

p(r,z,t )  =  J  p(r>z,f)e~l2n/tdf 
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Note  that  the  overall  phase  factor  e 


ilnf— 

c 

e  °  can  be  neglected  if  we  consider  that  the 


arrival  times  are  given  in  terms  of  “reduced  time”  T 


p(r,  z,  T)  = 


If 


vi f(r,z,f)e~i2KfTdf. 


(15) 


This  means  that  the  time  domain  is  heterodyned  around  t  =  — ,  and  the  use  of  a  reduced 


time  does  not  have  any  influence  on  the  autocorrelation  function.  The  model  also 


heterodynes  the  signal  by  shifting  the  center  frequency  to  zero  (d.c.),  in  order  to  reduce  the 
computational  load  associated  with  large  transform  sizes.  This  procedure  does  not 
introduce  any  consequence  to  our  algorithm. 

The  MMPE  model,  a  FORTRAN®  code,  is  the  latest  version  of  what  was  formerly 
known  as  the  University  of  Miami  Parabolic  Equation  (UMPE)  model  (Smith  and  Tappert, 
1993).  It  requires  a  series  of  input  files  with  strict  formats,  which  properly  represent: 

•  environmental  data  (sound-speed  profile,  bathymetry  of  water/bottom  interface, 
bathymetry  of  the  deep  layer  beneath  the  water/sediment  interface,  acoustic 
parameters  of  this  “deep”  layer,  and  acoustic  parameters  of  the  medium  just  below 
the  water/bottom  interface); 

•  source  data,  with  two  types  available,  wide-angle  and  vertical  line  array  (steering 
allowed).  We  only  used  the  wide-angle,  which  approximates  the  point  source.  The 
center  frequency,  bandwidth,  and  number  of  discrete  frequencies  (a  power  of  two, 
because  of  the  FFT)  are  also  required. 

All  these  input  files  are  specified  in  a  main  one,  the  “pefiles.inp”,  where  several  other  data 
(name  of  output  files,  grid  size  and  grid  steps,  c0,  size  of  vertical  FFT,  etc.)  are  contained. 
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The  output  is  presented  in  a  single  binary  file  composed  of  a  header  (with  the  information 
needed  for  post-processing),  and  the  complex  PE  field  function  \|/  at  selected  grid  points 
for  every  discrete  frequency. 
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III.  THE  FREQUENCY-DOMAIN  AUTO  CORRECTION  MATCHING 
(FACM)  LOCALIZATION  ALGORITHM 


The  algorithms  studied  in  the  Transient  Localization  Project  at  Naval  Postgraduate 
School  are  designed  for  data  from  a  single  receiver  hydrophone  and  a  point  source.  The 
basic  concepts  of  the  correlation  functions  are  outlined  in  Bendat  and  Piersol  (1971),  and 
the  routines  are  described  in  Miller  et  al.  (1996),  Pierce  (1996),  and  Hager  (1997).  From 
Brune  (1998)  we  have  that: 

Localization  algorithms  may  be  considered  generalized 
beamformers  in  which  the  plane  wave  replicas  have  been  replaced  by 
more  complicated  replicas  of  the  acoustic  propagation  (e.g.,  modes, 
beams,  or  the  vertical  pressure  field).  The  algorithms,  usually  referred  to 
as  processors,  are  in  most  cases  based  on  a  Hermitian  quadratic  product. 

The  exact  form  is  determined  by  the  constraints  that  are  put  on  the 
processor  output. 

As  discussed  in  Chapter  n,  replicas  are  simulated  received  signals  of  synthetic  sources 
located  at  the  points  of  a  search  grid.  Generated  by  a  propagation  model,  they  are  matched 
with  the  received  signal.  Considering  that  the  source  functions  are  equal  and  a  perfect 
model  is  used,  the  exact  match  should  occur  where  real  and  synthetic  source  .locations 
coincide.  An  ambiguity  surface  usually  represents  the  localization  algorithm  results. 

The  reciprocity  principle  states  that,  in  a  time-invariant  and  homogeneous 
environment,  the  acoustic  pressure  at  location  A  due  to  an  omnidirectional  source  located 
in  B  is  equivalent  to  that  of  location  B  due  to  an  omnidirectional  source  located  in  A. 
Using  this  principle,  we  reduce  the  computational  load  required  to  obtain  localization. 
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since  we  are  able  to  set  the  propagation  model  to  generate  replicas  from  the  receiver 
position  to  all  source  possible  locations  in  the  grid. 

In  this  chapter,  we  discuss  the  derivation  of  the  frequency-domain  autocorrelation 
matching  (FACM)  localization  algorithm.  Autocorrelation  matching  was  originally 
designed  for  use  in  the  time-domain  (TACM),  therefore  it  is  more  convenient  to  start  its 
study  treating  the  problem  in  time,  and  then  to  convert  the  results  to  the  frequency  domain. 
Previous  investigations  done  by  de  Kooter  (1997)  showed  that  the  TACM  algorithms 
present  very  small  footprints  and  large  phase  errors  at  higher  frequencies,  being  useful  for 
very  low  frequency  signals.  It  was  also  observed  that  matching  the  time-domain 
amplitude-squared  signals,  which  seemed  relatively  similar,  could  be  more  effective  and 
robust.  Therefore,  as  we  will  describe  in  the  following  paragraphs,  an  almost  natural 
further  development  led  to  the  frequency-domain  algorithm.  We  begin  with  a  complex 
pressure  time  series  P(t),  which  represents  a  detected  transient  arrival. 

The  signal  autocorrelation  in  the  time-domain  can  be  expressed  by 

Tpp{%)  =  J  P*(t)P(t  +  x)dt  (16) 

— o© 

or,  in  terms  of  the  frequency-domain  response,  by 

Tpp{x)  =  f  P*(f)P(f)e~i2nfxdf,  (17) 

— oo 

where 

P(f)  =  f  P(t)ei2%ftdt.  (18) 

— oo 
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Then,  the  normalized  autocorrelation  becomes 


Cpp(x)  — 


Tppi^) 

fWV' 


A  similar  derivation  applied  to  the  predicted  replica  signal  R(t)  leads  to 


(19) 


Crr(i) 


] [rr (t) 

f  I  *dV 


(20) 


The  autocorrelation  matching  algorithm  is  based  upon  the  inner  product  of  these  two 
quantities  and,  normalized,  in  time-domain,  it  can  be  expressed  by 


o© 

J  CRR(x)CpP(x)dx 


Arp  - 


f  \CRR^)\2dx\  \CPP(x)\2dx 

—  — oo  — oo  - 


1/2* 


(21) 


By  analogy,  in  the  frequency-domain  it  becomes 


f  CRRW)CppW)d$ 

Arp  =  — - — - — - 

>r  \cpp($)\2d<$> 

— oo 

where  Cyp  and  CRR  are  the  autocorrelations  in  frequency-domain  of  the  transient  signal 
and  of  the  replica,  respectively. 


Y  *  (22) 
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Since 


J  C"(<|))C~~(<^<t>=  J  <*|>[)  dt\R(tfeiln<i{]^\  dt'\P(t'fei2n¥ 


— oo  I— 


=  J~  dt'\P(t')\2f  dtmtfff  d§elimt'~t)~\  (23) 

— oo  — oo  1—  — oo  -J 

=  j~  \R(t)\2\P(t)\2dt 


Eq.  22  can  be  rewritten  as 


f  \R(tf\P(tfdt 


Arp  = 


f  \R(.t)\*dtf  \P{t)\Adt 


1/2' 


(24) 


In  order  to  reduce  the  influence  of  noise,  we  shall  remove  the  zero-lag  component  of 
the  autocorrelation  function,  which  is  equivalent  to  removing  the  mean  of  the  squared 
amplitude  of  the  time-domain  signal.  Therefore,  we  define 


|^(t)|2  =  |*(0|2-<l*W|2>r 


and 


(25) 


\P\tf  =  \P(tf-(\P(tf)t, 
where  <  >t  is  an  average  over  all  times.  Eq.  24  then  becomes 


(26) 


f  \R\tf\P\t)\2dt 


Arp  = 


il/2 


\f  \R'(t)\*dtf  \P'(t)\Adt 1 

L  —oo  — OO  J 

which  is  focused  on  the  influence  of  matching  the  non-zero  lag  components. 


(27) 
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However,  it  still  presents  an  ambiguity  in  (absolute)  time  which  was  not  present  in 
TACM,  and  we  shall  introduce  a  search  parameter  x  to  allow  us  to  slide  the  replica  signal 
in  time,  searching  for  the  optimal  match.  The  FACM  function  then  becomes 


f  \R’(t+ztf\P'vfdt 


j"  \R'(n\'dtf 

-  — oo  — 


(28) 


where  max%  is  the  maximum  value  of  the  function  for  any  value  of  the  search  parameter 
x. 


A  MATLAB  ®  implementation  of  the  algorithm  based  on  Eq.  28  was  used  to  generate 
all  the  ambiguity  surfaces  shown  in  this  work.  The  sliding-x  operation  was  performed 

according  to  the  following  development.  Let  g(t)  =  \P'(t)\2  and  h(t)  =  |i?’(t)|2,  then 

g(0  =  f  e~i2nftG(f)df  ,  (29) 

— oo 

and 

h(t)  =  f  e~i2nftH(f)df 


h(t  +  x)  =  J  e~iZrcf{t+x)H(f)df=  J  e~i2nftH(f)e~i2nfzdf  .  (30) 

— oo  —  oo 
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Therefore,  the  numerator  of  Eq.  28  becomes 


A(t)  =  J  h(t  +  x)g(t)dt . 


(31) 


If  h(t)  is  a  real  function. 


A(x)  =  J  h*{t  +  x)g{t)dt 


-i 


J  ei2KftH*(f)ei2%fxdf 


J  e~i2%rtG{f')df 


dt 


II 

ff'l 


i'2tc/t  f  -i2n(f'  - f)t 


H*(f)G(f')ei2Kfz  J  e 


dt 


df'df 


=  J  H*(f)G(f)ei2K/xdf  ,  (32) 

— oo 

knowing  that 

OO 

Je-i2^-^  =  5(f_/).  (33) 

— oo 

Hence,  the  search  over  x  was  in  fact  computed  by  a  multiplication  in  the  frequency 
domain  and  an  inverse  FFT  operation  represented  by  the  numerator  A(x).  Prior 
investigations  (Brune,  1998  and  de  Kooter,  1997)  had  shown  that  the  FACM  algorithm  is 
better  suited  for  larger  bandwidths  and  higher  frequencies,  when  compared  to  its  analog  in 


22 


the  time  domain,  the  TACM  algorithm.  Also,  FACM  is  more  robust  against  environmental 
mismatch  and  presents  better  peak-to-sidelobe  levels  and  larger  footprints. 
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IV.  SIGNAL  PROCESSING  AND  ENVIRONMENTAL 

ASSESSMENTS 

A.  MATCHED-FILTER  PROCESSING 

Matched  filtering  (MF)  is  a  process  where  the  received  signal  r(t)  is  correlated  against 
a  replica  of  the  transmitted  signal  s(t),  in  order  to  obtain  a  better  signal-to-noise  ratio 
(SNR)  and  reduce  the  effective  temporal  structure  of  the  source  down  to  an  ideal,  coherent 
pulse.  It  is  also  known  as  a  conjugate  filter,  since  its  frequency  response  is  basically  the 
conjugate  of  the  (transmitted)  signal  spectrum  (Tolstoy,  A.,  1993).  In  this  work,  matched- 
filtering  was  used  first  to  identify  the  short  duration  bursts  (0.1s)  in  the  2-min.  long  data 
files.  Second,  but  more  importantly,  to  process  the  bursts  and  obtain  data  emulating  the 
existence  of  a  coherent  source,  instead  of  the  actual  non-coherent  one.  This  also  provides  a 
baseline  for  the  localization  results  since  the  model-generated  replicas  are  based  on  a 
coherent  source.  Therefore,  this  is  equivalent  to  localizing  a  transient  with  known  source 
function. 

The  process  was  implemented  digitally  in  MATLAB® ,  according  to  the  expression 

Pit)  =  f  S*(f)R(f)ei2nfidf,  (34) 

- c« 

where  p{t)  is  the  matched-filtered  output  in  time-domain,  S*(f)  is  the  complex  conjugate 
of  the  transmitted  signal  in  the  frequency  domain,  and  R(f)  is  the  received  signal  in  the 
frequency  domain. 
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An  example  of  the  MF  output  in  time  domain  can  be  visualized  in  Figure  4.1,  and  its 
frequency  spectrum  in  Figure  4.2.  The  raw  data  is  also  shown  and,  by  comparison,  we  can 
observe  the  significant  noise  reduction. 


MafcMd-llierK)  buw  -Hy4#22,  Swruno  9 


Tnufsac]  Tune  {s] 


Figure  4.1.  Examples  of  a  Raw  (Left)  and  an  MF  Burst  (Right)  in  the  Time-Domain. 


Fr^uency^z] 


Figure  4.2.  The  Frequency  Spectra  of  a  Raw  Burst  (Left)  and  of  an  MF  Burst  (Right). 


26 


B.  PLANE- WAVE  BEAMFORMING 


Plane-wave  beamforming  is  one  of  the  many  techniques  used  to  distinguish  particular 
features  of  a  received  signal  while  reducing  the  influence  of  noise  -  which  is  usually 
generated  near  the  surface,  and  tends  to  occur  in  high  modes  or  high  angles  of  propagation 
(Jensen  et  al.,  1994).  It  consists  of  decomposing  the  signal  into  plane  waves  by  summing 
the  outputs  of  spatially  distributed  sensors  -  in  this  work,  eight  hydrophones  for  each  of 
the  three  sub-arrays  in  the  VLA. 

In  underwater  acoustic  propagation,  much  of  the  information  about  the  environment 
and  the  relative  source/receiver  location  is  contained  in  the  vertical  structure  of  the 
acoustic  field.  Therefore,  the  ability  to  resolve  the  vertical  structure  of  the  arrival  at  the 
receiver  can  lead  to  estimations  of  the  corresponding  source  location. 

A  beamformer  is  especially  useful  because  it  allows  us  to  give  preference  to  one 
direction  of  propagation  over  another,  implementing  a  spatial  filter.  In  this  work,  it  was 
used  to  observe  the  dominant  beams  when  the  signal  was  arriving  at  the  VLA,  providing 
information  for  a  comparison  of  the  localization  algorithm  results  by  the  distinct  arrival 
angle  structures. 

An  excellent  discussion  about  the  beamforming  process  is  presented  by  Defatta  et  al., 
(1988)  and  its  summary  follows.  Consider  an  infinite,  homogeneous  medium  containing  a 
line  array  of  equally  spaced  elements  positioned  along  the  y-axis,  and  an  undetermined 
number  of  remote  sources.  The  output  of  an  element  located  at  the  origin  of  coordinates 
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due  to  the  Ith  source  is  defined  as  Sz(r) .  Assuming  plane-wave  propagation  in  the  far-field, 
a  source  from  direction  0;  outputs 


~  Si  t  + 


nrfsin(0z)\ 
c  y 


(35) 


where  n  is  the  array  element  index,  d  is  the  spacing  of  the  elements,  0^  is  the  angle  of 


arrival  from  the  Ith  source,  and  c  is  the  speed  of  propagation.  To  accomplish  beamforming. 


we  must  apply  weights  0^  and  time  delays  to  the  individual  element  signals,  and  then 
coherently  sum  them.  Therefore,  the  output  of  the  beamformer  in  the  direction  0^  is 


N-  1 

bm{t)  =  X  wz5z|f  + 


nd[sin(0z)-  sin(0OT)] 


(36) 


n  =  0 


To  derive  the  directional  response  characteristics  of  the  beamformer,  we  shall  consider  a 
continuous,  complex  plane-wave  signal  propagating  across  the  array, 


ICO,? 


Si  -  e 


(37) 


Due  to  the  arrival  angle  0^,  the  sensor  outputs 


i(0,r  inklds'm(Q[) 

en(t)  =  e  e 


(38) 


2rc 

where  k,  =  t—  =  —  is  the  wave  number  of  the  Ith  source. 

1  h  c 
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The  beamformer  output  is  then 


where 


N-  1 


i(Ott 

e 


x 


wne 


m^/rf[sin(0/)-  sin(0m)] 


=  W(Q,)e 


i  (£>,!■ 


n  =  0 


(39) 


At-  1 

X 

n  =  0 


-  sin(0m)] 


(40) 


is  the  spatial  Fourier  transform  of  the  array  weights. 

The  beamforming  operation  applied  in  this  work  was  implemented  using  FFT  in 


MATLAB  ® .  Single  frequencies  considered,  the  time-delay  operations  described  above 
are  equivalent  to  phase  shifts  that  are  linear  functions  of  the  element  index.  FFT 
beamforming  is  the  implementation  of  these  phase  shifts  through  FFT  operations.  The 
processed  results  for  all  three  sub-arrays,  and  the  three  distinct  scenarios  can  be  seen  in 
Figures  4.3,  4.4,  and  4.5.  They  describe  the  physical  angle  of  arrival  as  a  function  of  the 
reduced  time,  and  provide  an  easy  visualization  of  the  multipath  nature  of  the  burst 
propagation.  The  results  for  sub-array  2,  however,  are  noticeably  affected  by  the  lack  of 
data  from  hydrophones  17  and  18.  We  can  verify  in  Figure  4.3  (scenario  11,  range  -500 
m)  that  the  energy  is  not  as  concentrated  in  the  lower  angles  as  in  the  other  two  figures 
corresponding  to  larger  ranges.  However,  in  all  of  them  we  observe  the  early  arrivals  of 
the  lower  angles,  the  short  period  in  which  the  energy  is  contained,  and  the  distinct 
multipath  structure. 


29 


Figure  4.3.  Beamformer  Outputs  for  Scenario  1 1  (Range  460  m.  Depth  45.7  m). 
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(c)  Sub-array  1 

Figure  4.5.  Beamformer  Outputs  for  Scenario  9  (Range  4650  m.  Depth  45.7  m). 
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C.  MODAL  ANALYSIS 


Normal  modes  are  another  physical,  frequency-dependent  representation  of  the 
underwater  acoustic  wavefield.  For  a  given  sound  channel,  the  acoustic  modes  represent  a 
unique,  natural  set  of  standing  waves  oscillating  along  the  depth  axis  (Chiu  and  Ehret, 
1994). 

The  method  involves  solving  a  depth-dependent  equation,  and  the  complete  acoustic 
field  is  then  constructed  by  summing  up  contributions  of  each  of  the  modes  weighted  in 
accordance  to  the  source  depth.The  modal  equation 


r  l  wjzh 

r  2  i 

0)  2 

;Lp(z)  dz  J  + 

2  Krm 

lc  (z)  J 

VJz)=  0, 


(41) 


with  boundary  conditions  VF(0)=  0  (pressure-release  surface  at  z=0) ,  and  —  =  0 

dz  Z  =  D 

(perfectly  rigid  bottom  at  z  =  D)  is  a  classical  Sturm-Liouville  eigenvalue  problem,  whose 

properties  are  well  known  (if  we  assume  that  p(z)  and  c(z)  are  real).  Some  are: 

•  the  modal  equation  has  an  infinite  number  of  solutions  which  are  like  the  modes  of 
a  vibrating  string; 

•  the  modes  are  characterized  by  a  mode  shape  function  ^m(z)  (an 
eigenfunction)  and  a  horizontal  propagation  constant  krm  (eigenvalues).  These 
constants  are  all  distinct  and  analogous  to  a  frequency  of  vibration;  and 

•  the  modes  are  unique  (orthogonal)  and  form  a  complete  set,  which  means  that  we 
can  represent  an  arbitrary  function  as  a  sum  of  the  normal  modes. 
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Based  upon  these  and  through  a  modal  decomposition  implemented  in  MATLAB,  we 
were  able  to  quantify  the  amplitudes  of  the  first  30  modes  as  functions  of  the  signal 
(discrete)  frequencies  and,  by  EFFT,  of  the  reduced  time.  The  SSP’s  used  were  the  ones 
discussed  in  Chapter  I  and  depicted  in  Figure  1.6,  applied  to  their  respective  cases.  A 
visualization  of  the  first  ten  modes  of  the  assumed  waveguide  for  a  frequency  f  =  900  Hz 
is  presented  in  Figure  4.6.  The  final  results  of  the  modal  decomposition  of  the  received 
signal  for  each  studied  scenario  are  shown  in  Figures  4.7,  4.8  and  4.9.  Note  the  energy 
concentration  in  the  lower  modes,  especially  scenarios  9  (5000  m)  and  16  (2500  m).  This 
is  consistent  with  stripping  of  higher  modes  at  longer  ranges.  The  energy  observed  in  the 
higher  modes  must  be  treated  with  skepticism  since  the  array  does  not  sample  those  modes 
very  well.  The  lack  of  tilt  information  in  this  analysis  also  makes  the  higher  mode 
decomposition  suspect. 
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Figure  4.6.  The  First  10  Normal  Modes  (f  =  900  Hz).  The  Bottom  Depth  is  89.9  m. 
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Figure  4.7.  Modal  Amplitudes  vs.  Reduced  Time  -  Seen.  11  (Range  460  m,  Depth  45.7  m). 


0.4  0.45  0.5  0.55  0.6  0.65  0.7 

Reduced  Time  [sec] 


Figure  4.8.  Modal  Amplitudes  vs.  Reduced  Time-Seen.  16  (Range  2194  m,  Depth  27.4  m). 
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Figure  4.9.  Modal  Amplitudes  vs.  Reduced  Time-Scen.9  (Range  4650  m,  Depth  45.7  m). 
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D.  SPATIAL  COHERENCE  (CORRELATION) 

In  order  to  analyze  the  spatial  coherence  of  the  signal’s  vertical  structure  at  the  VLA, 
single  records  from  each  of  the  three  scenarios  were  considered.  We  computed  the  peak 
correlation  of  each  phone’s  data  with  a  reference  phone  (hyd  #1,  closer  to  the  bottom;  16, 
mid-array;  and  32,  closer  to  the  surface)  normalized  such  that  the  autocorrelation  of  the 
reference  phone  data  was  unity.  The  results  provided  information  about  how  unique  the 
signal  was  across  the  space  (vertical  direction).  We  can  observe  in  Figure  4.12,  from 
scenario  9  (range  -5000  m),  that  the  signal  correlation  decreases  to  values  around  0.5  over 
about  5  meters  in  depth  for  the  reference  hydrophones  1  (bottom)  and  32  (top),  which 
corresponds  to  one  element  spacing  and  approximately  two  to  four  wavelengths.  The  same 
effect  can  be  verified  in  Figure  4.11  which  displays  data  from  scenario  16  (range  -2500 
m).  These  results  suggest  that  distinct  multipath  arrivals  at  any  given  time  are  interacting 
at  different  locations  over  the  array  length.  For  scenario  1 1  (range  -500  m)  runs,  however, 
the  cross-correlation  is  not  as  low  for  the  same  difference  in  depth,  as  shown  in  Figure 
4.10.  This  is  a  consequence  of  the  short  distance  between  source  and  receiver  which  does 
not  allow  large  variations  in  the  signal  nor  many  multipath  arrivals.  Also,  by  observing 
Figure  4.12,  we  notice  that  hydrophone  16’s  signals  are  more  correlated  than  the  ones 
,  from  hydrophones  1  and  32  for  scenario  9.  This  could  be  a  consequence  of  its  relative 
position  (mid-array),  but  the  same  effect  was  not  verified  in  scenarios  11  and  16. 
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Scenario 


(a)  Reference  Hydrophone  #  32  (Depth  28  m) 

Phone- to-phone  peak  oorreial ion  lor  Scenario  11  -MF 


(b)  Reference  Hydrophone  #16  (Depth  54.8  m) 


(c)  Reference  Hydrophone  #  1  (Depth  82  m) 

Figure  4.10.  Phone-to-Phone  Peak  Correlation  -  Scenario  11  (Range  460  m). 
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Phone-to-pbone  peak  correlation  lor  Scenario  16  -  MF 


(a)  Reference  Hydrophone  #  32  (Depth  28  m) 


(b)  Reference  Hydrophone  #16  (Depth  54.8  m) 


(c)  Reference  Hydrophone  #  1  (Depth  82  m) 

Figure  4.11.  Phone-to-Phone  Peak  Correlation  -  Scenario  16  (Range  2194  m). 
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Phone-lo-phone  peak  correlation  tor  Scenario  9  -  MF 
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(a)  Reference  Hydrophone  #  32  (Depth  28  m) 


Phone-lo— phone  peak  correlation  tor  Scenario  9  -  Mf 
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(b)  Reference  Hydrophone  #16  (Depth  54.8  m) 


Phone-to-pbone  peak  correlation  tor  Scenario  9  -  MF 
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(c)  Reference  Hydrophone  #  1  (Depth  82  m) 

Figure  4.12.  Phone-to-Phone  Peak  Correlation  -  Scenario  9  (Range  4960  m) 


V.  LOCALIZATION  RESULTS 


This  chapter  presents  a  compilation  of  results  obtained  after  a  series  of  runs  of  the 
FACM  algorithm.  For  every  scenario,  distinct  receiver  depths  were  tested  by  selecting  the 
hydrophone  data  to  be  used. 

Regarding  the  environmental  parameters  used  as  inputs  to  the  MMPE  model,  only  the 
sound  speed  profile  varied  according  to  the  respective  scenario.  We  used  only  the 
bathymetry  along  the  North  radial,  which  is  an  approximate  89.9  m  isobath  with  bottom 
composed  mainly  of  medium  to  coarse  grained  sands.  This  type  of  sand  has  the  following 

properties:  compressional  sound  speed  of  -1800  m/s,  sound  speed  gradient  of  25s"1,  shear 

speed  of  -250  m/s,  density  -2.03  g/cm3,  and  compressional  and  shear  wave  attenuation 
coefficients  about  0.1  dB/m  and  1.0  dB/m,  respectively  (Smith  et  al.,  1998).  During  the 
experiment,  the  data  was  recorded  at  a  sampling  frequency  of  4934.0  Hz.  For  processing, 
we  extracted  the  received  hydrophone  response  to  the  0.1s  bursts  in  files  of  4210  points 
which  corresponded  to  0.8533  seconds.  The  reason  for  this  procedure  is  that,  in  the 
MMPE  model  output,  we  have  512  discrete  frequencies  over  a  600  Hz  bandwidth. 

Therefore,  A/  =  1.171875  Hz,  and  we  match  T  =  ^  =  0.8533.  After,  we  processed  it 

again  to  reduce  its  length  to  2048  spectral  points  over  a  bandwidth  of  600  Hz  to  1200  Hz. 
The  model  output  was  zero-padded  (both  ends)  to  also  contain  2048  points. 
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The  spatial  grid  was,  in  most  of  the  runs,  128  points  in  depth  vs.  300  points  in  range.  In 
general,  it  corresponded  to  a  resolution  of  0.78  m  in  depth  and  10  m  in  range.  For  all  cases, 
the  source/receiver  relative  position  proved  to  be  a  very  important  factor  to  the  quality  of 
the  results.  This  can  be  observed  in  Figures  5. 1  (a)  and  (b),  where  the  ambiguity  surfaces 
for  scenario  9,  based  on  hydrophone  16  (a),  and  hydrophone  2  (b),  are  shown.  The  source 
is  located  at  depth  45.7  m,  range  4650  m,  and  its  location  is  depicted  as  a  white  dot.  While 
the  algorithm  with  hydrophone  16  (depth  =  54.8  m)  data  was  able  to  localize  it  within  the 
some  tolerance,  the  run  with  hydrophone  2  (depth  =  76.3  m)  data  did  not  provide 
localization  at  all. 

By  processing  the  matched  filtered  signal,  we  simulate  the  (idealistic)  situation  of  a 
coherent  source,  and  it  improved  the  localization  results  a  little.  However,  even  when 
processing  the  raw  data,  we  obtained  satisfactory  results  where  the  main  differences  were 
the  bigger  footprint  and  the  smaller  values  of  the  ambiguity  surface.  These  differences  can 
be  observed  in  Figures  5.2, 5.3,  and  5.4,  which  display  results  for  scenario  9,  scenario  16, 
and  scenario  11,  respectively. 

Observing  Figures  5.3(b)  and  5.4(b),  where  the  ambiguity  surfaces  from  runs  with  MF 
data  from  scenarios  16  and  11  are  shown,  respectively,  we  verify  that  the  FACM 
localization  algorithm  performed  reasonably  well  (errors  <  11%  in  range,  and  <  40%  in 
depth).  It  is  important  to  note  that  the  sound  speed  profiles  used  were  gathered  four  hours 
apart  from  experiments  of  scenario  9,  more  than  two  hours  apart  from  scenario  ll’s,  and 
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more  than  10  hours  apart  from  scenario  16’s.  The  mismatch  between  the  actual  and  the 
modeled  environmental  parameters  contained  in  the  SSP  is  likely  significant. 

We  also  observed  that  the  localization  error  decreased  with  increasing  range.  This  was 
expected  since  the  autocorrelation  techniques  are  strongly  dependent  on  the  uniqueness  of 
the  multipath  structure.  The  error  obtained  for  scenario  9  (range  -5000  m)  runs  were  in  the 
order  of  2%  in  range  and  30%  in  depth.  However,  scenario  11  (range  ~500  m)  runs 
presented  errors  as  large  as  20%  in  range  and  40%  in  depth. 
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(b) 

Figure  5.1.  Ambiguity  Surfaces  for  Scenario  9-(a)  with  Raw  Data  from  Hyd.16  (Depth 
54.8  m),  and  (b)  with  Raw  Data  from  Hyd.2  (Depth  76.3  m).  The  White  Dot  Depicts  the 

Actual  Source  Location. 
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Figure  5.3.  Ambiguity  Surfaces  for  Scenario  16-  (a)  with  Raw  Data  and  (b)  with  MF  Data, 
both  from  Hyd.  16  (Depth  54.8  m).The  White  Dot  Depicts  the  Actual  Source  Location. 
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Figure  5.4.  Ambiguity  Surfaces  for  Scenario  1 1-  (a)  with  Raw  Data  and  (b)  with  MF  Data, 
both  from  Hyd.  16  (Depth  54.8  m).The  White  Dot  Depicts  the  Actual  Source  Location. 
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VI.  CONCLUSIONS 


A.  SUMMARY  OF  FINDINGS 

Several  techniques  helped  us  to  assess  the  properties  of  the  acoustic  field  in  which  we 
were  applying  the  localization  algorithm  and  infer  their  influences  on  it.  Through  spatial 
coherence  analysis,  we  were  able  to  verify  that  more  unique  signals  led  to  better 
localization  results,  therefore,  ratifying  Brune  (1998).  Analyzing  the  beamformer  outputs, 
we  can  affirm  that  more  stratified,  discrete  structures,  corresponding  to  the  larger 
distances  between  the  source  and  receiver  also  provided  better  results.  Observing  the 
sound  speed  profiles,  we  identified  a  sound  channel  spanning  from  a  depth  of  about  30  m 
to  the  bottom.  With  the  modal  decomposition,  we  verified  that  the  lower  modes,  which 
were  closer  to  the  channel  axis  (~  45  m),  carried  more  energy.  When  selecting  a  receiver 
above  this  channel,  we  observed  poor  results  without  any  match,  and  the  best  results 
corresponded  to  hydrophones  positioned  closer  to  the  sound  channel  axis.  Another 
important  fact  to  note  is  that  the  depths  where  we  obtained  localizations  corresponded  to 
the  ones  containing  larger  amplitude  modes,  mainly  modes  two  to  six  depending  on  the 
source  frequency  chosen  to  be  sampled. 

Regarding  the  use  of  matched-filtered  or  raw  data  as  inputs  to  the  FACM  algorithm, 
our  work  suggests  that  it  is  not  a  major  issue  -  one  can  have  localization  results  applying 
raw  data  as  useful  as  the  ones  from  MF  data.  However,  these  results  were  based  solely  on 
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LFM  chirp  transmissions  and  may  not  apply  to  all  signals  generally.  Further  studies  are 
needed  to  examine  this  issue  more  thoroughly. 

B.  RECOMMENDATIONS  FOR  FUTURE  WORK 

Once  we  have  data  from  30  of  the  32  hydrophones  in  the  VLA,  an  almost  natural 
consideration  would  be  the  multiple-phone  analysis  of  the  localization  algorithms.  While 
computational  load  would  certainly  be  an  issue,  the  comparison  with  single-phone 
analysis,  and  possibly  better  results,  would  be  worth  it.  In  this  case,  one  should  consider 
the  use  of  the  array  tilt  information  contained  in  the  data  files  to  correct  the  arrival  angle, 
vertical  structure,  hydrophone  depth,  etc. 

Another  possible  approach  would  be  the  use  of  non-traditional  processes,  such  as 
artificial  neural  networks  (ANN)  and  minimum- variance  spectrum  estimations  (MVSE), 

.  both  in  the  arrival  structure  analysis  and  in  the  data  preprocessing.  This  could  surpass  the 
over-simplifications  of  the  traditionally  used  techniques  and  lead  to  better  results  (Ma, 
Y.L.,  1997).  Also,  since  the  shallow  water  environment  is  generally  a  dispersive  channel, 
and  presents  effects  such  as  strong  bottom  reverberation,  signal  distortion  and  noise 
fluctuation,  an  analysis  of  the  signal  mismatch  influences  could  introduce  improvements 
,  to  the  process. 


50 


Throughout  this  work,  only  three  scenarios  were  considered.  As  the  NUWC-TSP 
experiment  provided  at  least  14  other  sets  of  data  where  parameters  such  as  source 
location,  signal  bandwidth,  and/or  pulse  duration  differed  from  the  ones  we  studied,  it 
could  be  useful  to  process  them  to  try  to  confirm  our  results.  Moreover,  it  would  allow  a 
comparison  of  the  FACM  performance  for  longer  pulses  and  higher  frequencies. 
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